function [val] = myacf(x)
   
   alpha = 1.96;
   
   rho = @(xy,h) mygamma(xy,h)./mygamma(xy,0);
   n = length(x);
   nlag = floor(min(10*log10(n),n-1));
   
   vall = arrayfun(@(h) rho(x,h),1:nlag);
   x2 = x.^2;
%    mm = mygamma(x2,1);
%    C = mycov2(x-mean(x),1,1);
   var = 1+ arrayfun(@(h)mygamma(x2,h),1:nlag)./(mygamma(x,0).^2);
   band = -sqrt(var/n);
%    minval = 1.2*min(min(vall), min(-alpha*band),-alpha./sqrt(n));
%    maxval = 1.2*max(max(vall)',  max(alpha*band),alpha./sqrt(n));
   
   autocorr(x,nlag);
   hold all;
   plot(1:nlag,-alpha*band,'r');
   plot(1:nlag,alpha*band,'r');
   a = axis;
   axis([a(1:2) -0.5 1]);
   hold off;
   
end

function gamma = mygamma(x,h)    
    n = length(x);
    h = abs(h);       
    x = x - mean(x);
    gamma = sum( x(1:n-h).*x(h+1:end) )./n;    
end

function xx = mycov2(x,h,k)
% calculate the covariance(i,j): (x_{t})^2 .* x_{t-k} .* x_{t-h}
    
%     x = x - mean(x);

    n = length(x);
    x1 = (x(max(h,k)+1 : n).^2);
    x2 = x(1 : n-max(h,k) );
    x3 = x( max(h,k) + 1 - min(h,k) : n - min(h,k));
    
    xx = sum(x1.*x2.*x3)./n;
   
end